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ABSTRACT 

Galactic disks consist of both stars and gas. The gas is more dynamically 
responsive than the stars, and strongly nonlinear structures and velocities can 
develop in the interstellar medium (ISM) even while stellar surface density pertur- 
bations remain fractionally small. Yet, the stellar component still significantly 
influences the gas. We use two-dimensional numerical simulations to explore 
formation of bound condensations and turbulence generation in the gas of two- 
component galactic disks. We represent the stars with collisionless particles and 
follow their orbits using a particle-mesh method, and treat the g 9.1 S 9.1 S 911 isother- 
mal, unmagnetized fluid. The two components interact through a combined grav- 
itational potential, which accounts for the distinct vertical thickness of each disk. 
We ensure our results are not contaminated by particle Poisson noise, which can 
seriously compromise simulations of self-gravitating systems. Using stellar pa- 
rameters typical of mid-disk conditions, we find that models with gaseous Toomre 
parameter Q g < Q S)Cr it ~ 1.4 experience gravitational runaway and eventually 
form bound condensations. This Q 9iC nt value is nearly the same as previously 
found for razor-thin, gas-only models, indicating that the destabilizing effect of 
live stars is offsets the reduced self-gravity of thick disks. This result is also 
consistent with empirical studies showing that star formation is suppressed when 
Q g £ 1 — 2. The bound gaseous structures that form have mass 6 x 10 7 M Q each, 
~ 10 times larger than the thin-disk Jeans mass for a gas-only disk; these repre- 
sent superclouds that would subsequently fragment into GMCs. Self-gravity and 
sheared rotation also interact to drive turbulence in the gas when Q g £ Q 9)Cr it- 
This turbulence is anisotropic, with more power in sheared than compressive mo- 
tions. The gaseous velocity dispersion is ~ 0.6 times the thermal speed when 



-2- 



Q g ~ Qg.crit- This suggests that gravity is important in driving ISM turbulence 
in many spiral galaxies, since the low efficiency of star formation naturally leads 
to a state of marginal instability. 

Subject headings: galaxies: ISM — galaxies: structure — instability — ISM: 
kinematics and dynamics — solar neighborhood — stars: formation — stars: 
kinematics 



1. Introduction 



Large-scale galactic disk evolution is governed by the interaction between self-gravity 
and other reinforcing and/or countervailing forces. Some of these forces are associated with 
sheared galactocentric rotation, and others are associated with random velocities (thermal 
and turbulent) of the gas and stars. Under certain circumstances, gravitational instabilities 
may grow and radically alter the disk structure as they develop nonlinearly. This is believed 
to be how giant molecular clouds (GMCs), and subsequently star-forming H II regions, 
originate. Under other circumstances, the overall disk morphology is not changed, but self- 
gravitating modes grow and interact in such a way that significant energy is transfered from 
ordered to disordered forms. This "disk heating" may in turn affect the system's subsequent 
susceptibility to self-gravitating instabilities. 

Because the gas component is much cooler than that of the stars, and because gas can 
radiate energy as it nonlinearly compressed, the effects of self-gravitating instabilities are 
much more pronounced in the interstellar medium (ISM) than in the stellar disk. The stel- 
lar component can nevertheless be quite important to the initial growth of self-gravitating 
modes (or waves), because it contains so much of the disk's mass (~ 75 — 90%, compared to 
~ 10 — 25% in gas). While the stars are essentially collisionless, most previous work on grav- 
itational instability in two-component disks has treated both star s and gas as two i sotherm al 
fluids with differen t sou nd speeds. This approach was taken by |Jog fc Solomon Jl984al lbL 
Elmegreenl (I1995b[) . and iJogl (119961 ) in studying axisymmetric self-gravitating modes, and 
by Jogl ( 19921 ) in studying non- axisymmetric self-gravitating modes, for combined star /gas 
disks. Adopting a f ull kinetic treatme n t for the st ellar componen t in two-component disks, 
on the other hand, iLin. Yuan, fc Shul (119691 ) and iRafikovl (120011 ) used linear perturbation 
theory to derive dispersion relations for non-axisymmetric WKB waves and axisymmetric 
modes, respectively. These linear-theory analyses cannot capture, however, the ultimate 
outcomes - including GMC formation and turbulent driving - which result when localized, 
self-gravitating structures grow; quantifying this development requires nonlinear numerical 
simulations. 
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The main class of self-gravitating instability that develops under conditio ns of strong 



sheared rotation (in outer galaxies, away from spiral arms) is the swing amp lifier (IGoldreich fc Lynden-Bell 
1965bl ; lJulian fc Toomrd Il966l ; iToomrd Il98ll ). In iKim fe Ostrikerl (120011 ) (hereafter Paper 
I), we identified the basic requirements for swing to occur, and reported on the results of 
extensive numerical simulations of swing amplification in gaseous disks with or without mag- 
netic fields. Although in linear theory the amplification magnitude varies continuously with 
respect to the Toomre parameter Q g (see eq. [6] for definition), Paper I showed that non- 
linear interactions among swing-amplified density filaments or wavelets lead to a threshold 
phenomenon. Disks with Q g > Q 9>C rit remain stable, while disks with Q g < Qs.crit experience 
gravitational runaway, forming bound condensations. For the razor-thin disk models of Paper 
I, we found that Q 5)C rit ~ 1.2 — 1.4, with the largest value resulting from models with subther- 
mal magnetic fields. Given the similarity between the numerically-obtained Qn. ^t. values and 



observationally i nferred thresholds for active star formation in external galaxies (I Quirk! 11972 



Kennicuttlll989l : ICaldwell et aHll99ll : iMartin fe Kennicuttll200~ll : IWong & Blitzll2002l ). Paper 
I supported the notion that self-gravitating instabilities define the star-formation "edges" of 
disks. 

While successfully demonstrating the nonlinear threshold behavior of the swing ampli- 
fier, and quantitatively yielding good agreement with observations, the models considered in 
Paper I suffered from two important drawbacks. First, by assuming infinitesimally thin disks, 
they overestimated self-gravity, which tends to overestimate <3 giCr i t ; second, they did not al- 
low for a dynamically- active stellar disk, which tends to underestimate Q ffiCr it- Subsequent 
simulations of pure-gas disks in three dimensions confirmed the basic nonlinear threshold 
behavior, but also showed that for realistic disk temperatur es the nonzero thickness of th e 
disk can indeed significantly lower the value of Qg )C ?\t to < 1 (IKim. Ostriker. fc Stone! 120021 ). 



In this paper, we address the second limitation of Paper I by studying nonlinear evolution 
of gravitational instabilities in two- component galactic disks, consisting of both stars and gas. 
We distribute collisionless particles to represent stars and follow their orbits using a direct N- 
body method, while adopting a hydrodynamical approach for gas. The two components are 
allowed to interact with each other through the combined gravitational potential (calculated 
on a single mesh). We also allow for nonzero (and differing) vertical thickness of both stellar 
and gaseous disks in solving the Poisson equation for each component. 

Given that finite thickness effects tend to reduce Q s , cr it values below observed star for- 
mation thresholds in gas-only models, a compelling question is whether allowing for an active 
stellar disk compensates for this effect. Thus, one of our primary objectives is to quantify 
the critical Q g value for gravitational runaway when both a live stellar component and a 
thick disk are explicitly considered. A second important goal is to understand what happens 
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when Q g is sufficiently large that gravitational runaway does not occur. In particular, the 
large-scale flows that are driven when gravity is strong (but not too strong) may be a key to 
powering turbulence within the disk. 

Many processes may contribute to turbulence in the diffuse gas. In addition to self- 
gravitating instabilities, other candidate mechanisms that have been proposed include ther- 
mal instabilities, the Parker and magnetorotational instabilities, dynamical instabilities in 
spir al shocks, stellar winds, and expanding stellar H II regions and supernova blast waves 



; sec 



Mac Low fc Klessenl I2004J ; lElmegreen fc Scald 120041 for recent reviews) . Traditionally 



stellar processes have been presumed to be the dominant turbulence drivers, but the lack 
of observed spatial c orrelation between yelocity dispersions of H I gas and regions of high 
mass star formation (jDickey et al.lll990l ; Ivan Zee fc Brvantlll999l) a rgues that other s ources 
must be important as well. In particular, lElmegreen! (120021 ) and lElmegreen et al.l (120031 ) 
have suggested that the swing amplifier generates large-scale motions that cascade to pro- 
duce turbulence and structure on small scal es as well as large scales. Initia l numerical 
studies of this process have been pursued by IWada. Meurer. fc Norman! (120021 ) . who per- 
formed two-dimensional, thin-disk galactic simulations with ISM cooling. They found that 
a clou dy medium develops, w i th a y elocity dispersion a v ~ 2 — 5 km s _1 . On the other 
hand, iKim. Ostriker. fc Stond (120031 ) found via three-dimensional isothermal simulations, 



with sound speed c g = 7 km s" 1 and Q g ^ 1.7, that turbulence was driven to levels less 
than 1 km s . Based on just these few models, it has not yet been clear how important 
"swing-driven" turbulence is more generally. 

In the present work, we determine how the turbulent velocity amplitude of the gas, a v , 
varies as a function of Q g , and also assess how the presence of a live stellar component affects 
a v . By running a number of models with a range of particle numbers, and also performing 
control models with "passively-evolving" stars, we are able to subtract out the contribution 
to a v due to Poisson noise in the particle distribution. We show that quite high particle 
resolution is in fact necessary for accurate determination of a v in the case when Q g is near 
Qg,crit] this is of practical importance because real disks naturally evolve toward near-critical 
states. 

This paper is organized as follows: In §2, we present the basic equations we solve for gas 
and stars, introduce "thick-disk" gravity, and describe our model parameters. In §3, we re- 
visit the linear theory for axisymmetric gravitational instability in two-component disks with 
an emphasis on the stabilizing effect of nonzero disk thickness (omitted in previous studies). 
In §4, we describe numerical methods we employ for our simulations, and present results of 
code tests. We demonstrate that large Poisson noise in the stellar particle distribution can 
lead to spurious results for iV-body systems that are nearly or marginally gravitationally 
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unstable. In §5, we report on our simulations of two-component disks. We present numeri- 
cal results on Q g thresholds for (non-axisymmetric) gravitational runaway, and quantify the 
overall level and spectral properties of gravity-driven turbulence when Q g £ Q 9iCr i t Finally, 
we summarize our results and discuss their astronomical implications in §6. 



2. Equations and Model Parameters 

2.1. Basic Equations and Thick-Disk Gravity 

In this paper we investigate the dynamics of local galactic disks composed of gas and 
stars. We follow the evolution of the gas by solving the hydrodynamic equations, while 
employing an iV-body method for the stars. The gas is assumed to be isothermal and un- 
magnetized, and we assume that all physical variables other than the gravitational potential 
are independent of the vertical height. As in Paper I, we consider a patch of the disk orbiting 
the galaxy with a fixed angular frequency fio = ^oz, and set up a local Cartesian frame with 
x and y referring to radial and azimuthal coordinates, respectively. We consider a simulation 
box with size Lx L. The equilibrium background velocity arising from galactic rotation rela- 
tive to the center of the box at x = y = is given by v = — g^o^y? where q = —d In Q/d In R 
denotes the dimensionless local shear rate. The basic dynamical equations we solve for the 
gaseous part in the local frame are identical to those presented in Paper I, except that we 
adopt here a "thick-disk" gravitational kernel in solving the Poisson equation (see below). 

To follow the evolution of the stellar component, we use collisionless particles. Restrict- 
ing particle motions to be in-plane, the shearing-sheet equations of motion relative to the 
center of the box are given by 

r = 2gn^x-2a xr-V($ 9 + $ s ), (1) 

where r = xSt + yy is the position vector of a particle, the dots den ote derivatives with 



respect to time (e.g.. lJulian fc Toomrelll966l ; Wisdom fc Tremaindll988l ). and $ 9 and $ s are 



the gravitational potentials from gas and stars, respectively. The spatial distribution of the 
particles at any given time will give the stellar surface density £ s , which in turn yields $ s 
via 

V 2 $ s = 47rG£ACz). (2) 

Here, h s (z) represents the vertical distribution of the stellar density and satisfies the normal- 
ization condition h s (z)dz = 1. The Poisson equation for the gas takes the same form as 
equation (j2J) with subscripts changed appropriately. 
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Studies of disk dynamics often assume h(z) = 5(z) for simplicity. This thin-disk ap- 
proximation would be valid as long as perturbations of interest do not critically rely on 
the vertical dimension and their wavelengths A are much larger than the disk scale height 
H. For waves with A approaching 2ttH, however, the thin-disk approximation overestimates 
self-gravity at the disk midplane. This can be particularly severe in a two-component sys- 
tem in which the characteristic length scales and disk scale heights for gas and stars are 
quite different from each other. For instance, the critical wavelengths for axisymmetric 
gravitational instability u nder solar neighborhood c onditions are ~ 1 kpc for stars (e.g., 
Binney fc Tremaind 119871 ) and ~ 1 kpc for gas (e.g.. lElmegreenlll995al ; Paper I). The stellar 
scale height H s ~ 330 pc (IChen et al.ll200ll ) is non-negligible compared with the wavelength 
of the most vulnerable modes in the combined system. Moreover, as perturbations grow via 
swing amplification they form thin overdense filaments with widths comparable to the gas 
disk's scale height (Paper I; see also §5.1|) . so that the thin-disk gravity may potentially lead 
to spurious fragmentation. It is thus desirable to allow for finite disk thickness in solving 
equation (j2J). 

Suppose that h s (z) does not vary with time and that E s is periodic in x and y (or linear 
combinations of coordinates; see §4.1). Using the Green's function method, one can show 
that equation (j2J) yields the solution at z = (omitting the subscript s) 



$(k) 



2vrGS(k) 



h(z')e- lkz ' l dz', 



(3) 



where $(k) and S(k) refer to the Fourier components with wavevector k. When h(z) 
e~w/ H /(2H), equation (JHJ) simply becomes 



$(k) 



2vrG£(k) 

|k|(i + |k|#y 



(4) 



which we refer to as "thick-disk" gravity, as opposed to thin-disk gravity for which H is set 
equal to zero. This result was obtained through a different technique by lElmegreenl (119871 ) 
and has been shown to be quite a good approximatio n for the reduction of self-grav ity in 
disks with various density profiles (e.g., iRomeol Il992l ; iKim. Ostriker. fc Stond l2002l ). By 
direct numeral integration, one can see that equation PJ slightly underestimates the rea l 
potential for purely self-gravitating disks, but not more than 14%FI IKim fc Ostriker ( 2006 ) 
further found that when equation (j3J) is used to calculate self-gravity, simulations of gas flow 



1 When equation ((4]) is used to obtain an axisymmetric dispersion relation for a disk with H g — c 2 g / (irGY} g ) 
(i.e. H g det ermined by gas self-gravity alone), the resulting Q cr it = 0.65 is very close to the value Q C rit = 0.68 
obtained bv lGoldreich fc Lvnden- Bell (|1965aj ) as an exact criterion for axisymmetric instability. 
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across two-dimensional spiral shocks produce results very similar to fully three-dimensional 
models. In what follows, we shall use equation (j3J) to compute the self-gravity of gas and 
stars with scale heights H g and H s , respectively. 



2.2. Model Parameters 



As initial models we consider a composite disk of gas and stars with uniform surface 
densities £ 9 o and £ s o- The gaseous disk is set to remain isothermal throughout the entire 
evolution with an effective speed of sound c g = 7 km s" 1 , corresponding to a mean Galactic 
thermal pr essure P/k ~ 2000 — 4000 K cm -3 (jHeilesll200ll ) and mean midplane density hh ~ 



0.6 cm 3 (IDickey fc Lockmanl 119901 ). similar to solar neighborhood values. For the stellar 
disk, we initially distribute particles according to the Schwarzschild distribution function 



f(vx,v y ) 



■°s0 



2tt<t s x u 



exp 



s.y 



21 



2cr 2 

w s,x 



2a 2 



(5) 



where a s>x and a SjV 



are the x- and y-components of the particle velocity dispersion, re- 
spectively. The values of a StX and a s>y will vary over time as the particles respond to the 
perturbed gravitational potential. Note that S s q = / fdv x dv y and cr Sjy /cr S:X = (1 — g/2) 1//2 in 



equilibrium. We initially adopt a S:X ■ 
neighborhood conditions with a flat 



30 km s 1 and a S}V = 2 1 / 2 cr. 9rT corresponding to solar 



1) rotation curve (jBinney fc Tremaindll987l ) 



When written in dimensionless form, the basic dynamical equations are characterized 
by the following parameters 



n.) 



K Q Cg 



go 



GTjgoL 



x 



3.36GS 



s0 



47r 2 GS 



(6) 



(7) 



.s0 



together with H g /L and H s /L. The Toomre Q parameters defined by equation ([6]) give the 
surface densities relative to the critical values at Q g = Q s = 1 for axisymmetric gravitational 
insta bility in a razor-thin, gas-only or star-only disk (IToomrei 1 19641 ; iBinney fc Tremaine 
19871 ). In equations (J7]), the rij (X) parameter is the ratio of the simulation box size L to the 
critical wavelength A 9iCr it = c 2 /GS 90 (A S)Cr it = 47r 2 GS s0 //to) f° r axisymmetric gravitational 
instability in a razor-thin, non-rotating gaseous (rotating, cold stellar) disk. 

Since the parameter space is very large and since we are primarily interested in the 
dynamical evolution of the gaseous component in combined disks, we vary only Q g (or 
equivalently S s0 ) and fix all the other disk properties. For the stellar part, we adopt solar 
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neighborhood values kq = 36 km s 1 kpc 1 (IBinnev &: Tremaine 1 19871). S,n = 35 pc 



(iKuiiken fc Gilmordll989h . and H s = 330 pc flChen et al.ll200ll ; Karaali et alJl2004r l. This 
gives Q s = 2.1, so that in our model disks the stellar component alone, even in the razor-thin 
limit, would be immune to axisym metric instability. We adopt H g = 170 pc for the effective 
scale height of the local ISM (e.g., iBoulares fc Coxl|l990l ). For the box size, we adopt X = 2 
or L = 9.2 kpc, which is large enough to resolve the most susceptible modes of the swing 
amplifier in stellar disks (lToomrelll98ll ); test simulations we have performed confirm that 
results are insensitive to X as long as X > 2. 

Our fiducial model has Q,„ = 1.4 corresponding to £ 9 o = 13 M pc -2 for the gas disk 



(e.g., iHolmberg &: Flynnl 120001 ). but we also allow for various values of Q g in the range 
of 1 - 3. These models will allow us to determine the critical Q g value for gravitational 
runaway and to study the characteristics and level of turbulence driven by self-gravity when 
Q g ^ Qc/,crit- Since equations (|SJ) and © give nj = 3.7Xa s , x / \c g Q g Q s ) ~ 15/ Q g for the 
adopted set of parameters, we see that our simulation box is large enough to contain at least 
5 Jeans wavelengths of the gaseous medium. We define the orbital period t or b = 2ti/Q = 
2.4 x 10 8 yr (Q /26 km s _1 kpc -1 ) -1 and use it as the time unit in our presentation. 



3. Axisymmetric Stability 



The gravita tional stability of two-c o mponent disk s to axisymmetric p ertur bations wa s 
first analyzed by I Jog fc Solomon! (11984al ) , iRomeol (119921 ) , lElmegreenl ( ]1995bl ) , and I Jog (119961 ) , 
who treated the stellar disk as an isothermal fluid. While a fluid description of stellar par- 
ticles is a good approximation for a range of wavenumbers in disks near the threshold, it 
generally fails when disks are further from instability or when the wa velengths of pertur- 
bations are sma ll com pared to the epicyclic excursions of stars (e.g., IBinnev &: Tremaine 



19871 ). iRafikovl (120011 ) instead used a collisionless description of the stellar population and 



derived a dispersion relation for axisymmetric waves in star-gas disks assuming that the disk 
is razor-thin. In this section, we revisit axisymmetric instability with a particular emphasis 
on the effect of nonzero disk thickness. The results of this section will be used to check our 
numerical technique in §4.21 



We refer the reader to IRafikovl (120011 ) for a detailed derivation leading to the thin-disk 
dispersion relation for two-component disks. We follow the same steps, except we use thick- 
disk equation for the potential and density pairs in place of Rafikov's equation (14). 
The resulting dispersion relation for axisymmetric modes in a two-component disk of finite 
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vertical extent reads 

2irGkZ g0 2nGkZ s0 F(u;/K, k 2 a 2 s jK 2 ) _ 

(n 2 + k 2 c 2 g -uj 2 )(l + kH g ) + (k 2 -uj 2 )(l + kH s ) ~ ' U 

where uj and k are the frequency and radial wavenumber of per turbations, respectively, an d 



T is the stellar "reduction factor" defined by equation (6-45) of iBinney fc Tremaind (119871 ) 



which accou nts for t he ef fects of the stellar velocity dispersion. A similar relation was 



presented in iRomeol (119921 ) who introduced a single effective scale height for the combined 
disk instead of using different heights H g and H s for each compo nent. In the l imit of 
H g = H s = 0, equation (jSj) recovers the thin-disk dispersion relation of iRafikov! (120011 ). Note 



that although a combined thick disk can be mapped to a razor-thin counterpart with lower 
surface density, the density reduction factors of the gaseous and stellar components differ 
from each other and cannot be determined a priori since they depend on the perturbation 
wavenumber. 

We determine the marginal stability condition by solving equation (jSJ) with uj 2 = 0, and 
plot the results in Figure [TJfor selected values of R = c g /cr S)X . Thic k curves a re for thick disks 



with scale heights H g = 170 pc and H s = 330 pc. Results from IRafikov! (120011 ) for razor- 
thin disks are also plotted as light curves for comparison. The parameter domain below 
each marginal curve corresponds to the stable regime. It is apparent that axisymmetric 
gravitational instability is mainly determined by the gaseous component for small R and 
Q s 1. When R is small, the characteristic unstable length scales in (standalone) gaseous 
disks are much smaller than those of the stellar disk. In this case, the relative dominance of 
the gaseous to stellar components changes rather abruptly as Q g and Q s vary along a given 
marginal stability curve, as manifested by the cusp in the R = 0.1 curve. As R increases, the 
distinction between the characteristic length scales in the gas and stars becomes smaller, and 
the marginal curve bends smoothly. A system with solar neighborhood parameters, marked 
by the dot at (Q s ,Q g ) = (2.1, 1.4), is highly stable to axisymmetric modes when treated as 
thick, while it would be marginally stable under the thin-disk approximation. 

For varying disk thicknesses, Figure [2] plots the Q 9jCr i t values, and the corresponding 
marginal wavelength of perturbations, A marg (A < A marg modes are stable). Self-gravity 
weakens when H g , H s , or Q s increases, tending to require lower Q gtCT it- A razor-thin gas-only 
disk has Q 9iCr it = 1.0. This value increases to Q 9)Cr it = 1-27 if the contribution from a thin 
stellar component with Q s = 2.1 is considered. For disks with thicknesses H g = 0.87c 9 /k and 
H s = 0A<j sx /k, with Q s = 2.1 (similar to the solar neighborhood conditions), however, we 
have Q g ,crit = 0.67 at A 

marg ~ 2.3 kpc, as marked by dots in Figure [21 This indicates that the 
stabilization caused by finite disk thickness is considerable, even larger than the destabilizing 
effect of stars. Solutions of equation (jSJ) show that when H g = 0.87c g /n and H s = 0Act SjX /k 
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are fixed, thick disks with Q g > 1.4 are stable to axisymmetric perturbations unless Q s < 0.8 
(which is a very unrealistic range). Therefore, the naive expectation that including a stellar 
component increases Q 9jCr it over unity fails in disks with realistic thickness. The fact that 
Qg, th ~ 1-4 for observed star formation thresholds in disk ga laxies has often b een attributed 
to gravitational instabilities in two-component systems (e.g.. lKennicuttlll997l ). However, our 
results show that axisymmetric modes cannot be providing the needed instability, when disk 
thickness is properly considered. We will show in §5 that instead, swing amplification of 
non- axisymmetric perturbations are in fact able to drive two-component thick disks with 
Q g £ 1.4 into eventual gravitational runaway. 



Numerical Method and Code Tests 



4.1. Numerical Method 



We have run a number of two-dimensional simulations for star plus gas systems, varying 
Q g and the number of stellar particles while fixing H g , H s , and Q s . All the models have 
256 x 256 resolution. We f ollow the nonlinear evol ution of the gaseous part using the same 
version of the ZEUS code (IStone fc Norman! 1 199 21 ) as in Paper I; the reader is refereed to 
Paper I for a detailed description of the code and the boundary conditions we adopt. In this 
subsection, we describe the numerical methods mainly for the stellar part. 

We use up to iVp tl = 2 x 10 6 particles to represent a stellar disk and evolve them based 
on a particle-mesh (PM) algorithm. With its high velocity dispersion, the stellar component 
has a characteristic Jeans scale much larger than that of gas, so the resolution limit imposed 
by the PM method is entirely tolerable. In addition, the relaxation time of the particles 
amounts to tn = erf x Ax/(iiG 2 T, s0 m) ~ 1 x 10 3 t or b for th e parameters w e adopt, where Ax 
is the grid size and m is the mass of individual particles (iRybickil 119711 ). This time is long 
enough to ensure that the evolution of the stell ar system in our PM code is not seriously 



influenced by particle noise and relaxation (e.g.. IWhitdll988l ). 



The initial distribution (eq. [5]) of the particle positions and velocities are realized 
using a pseudo-random number generator. At each time step, the stellar surface density 
£„ is calculated on a 256 x 256 mesh via the triangular-shaped-cloud assignment scheme 



( jHockney fc Eastwood! 119881 ) . We solve equation ([4]) by employing the fast Fourier trans- 
form method o n a sheared coordinate grid in which the surface density is perfectly periodic 
( lGammidl200ll ). The calculated potential is differentiated to give the gravitational force field 
at Cartesian mesh zone centers, which is then interpolated back to the particle locations. 
We integrate equation ([1]) using a modified predictor-corrector scheme which is second order 
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accurate in time (e.g.. lMonaghanlll989l ). To handle the kinematics of the background shear 



self-consistently, we apply the shearing-box boundary conditio ns in which ^-boundaries are 
perfe c tly periodic and the x-bo undaries are shearing-periodic (IHawley. Gammie. &: Balbus 
19951 ; iHuber fc Pfennigerl 1200 ll ). These boundary conditions are fully consistent with our 



local models; particles leaving the simulation box from one x-f&ce re-enter with shifted y 
positions though the opposite x-face. 



4.2. Code Tests 

We have checked our numerical methods on a number of test problems: one-dimensional 
or two-dimensional waves in gaseous or stellar disks separately or in combination. In this 
subsection we present the results of three tests that not only verify the accuracy of our 
numerical code, but also provide information on the necessary number of stellar particles to 
prevent numerical artifacts in the evolution of a self-gravitating stellar system. 

Figure [3] plots the test results of axisymmetric density waves in razor-thin, stellar disks. 
For these tests, we consider one-dimensional, star-only disks with H s — 0, q — 1, Q s — 1.1 
(close to the marginal state), and X = 1 or 2. For each run, N pt \ = 10 4 particles are used. 
Frequencies obtained from the spatial and temporal Fourier analyses of the surface densities 
are given as solid triangles and open circles for the X = 1 and 2 cases, respectively. The 
numerical results are within ~ 2.5% of the analytic predictions (solid lines) from equation 
([8]) for the chosen set of disk parameters. This confirms the stellar code's performance. 

The second test addresses axisymmetric gravitational instability of a two-component 
disk, for which we choose Q g = 1.2, Q s = 1.5, q = 1, H g = 0.05c 5 /k, H s = 0.1<j StX / k, 
and c g /a SjX = 0.4. The linear theory in §|3] predicts that such a combined disk is unstable 
with a maximum growth rate of ~ 0.894f2o, which occurs at A max = 3.53a S}X / n. Setting the 
one-dimensional box size L = A max and assuming the gas is adiabatic with index 7 = 5/3, 
we have run two models with differing iV pt i. Figure H] plots the resulting time histories of 
the maximum surface densities in both gaseous and stellar components. The dotted and 
solid lines correspond to N pt \ = 10 4 and 10 5 cases, respectively. Poisson noise in the stellar 
distribution immediately imposes (via gravity) strong density perturbations in the gaseous 
disk, which would otherwise be uniform. The perturbations in stars and gas soon become 
coherent and grow exponentially. The model with jV p ti = 10 4 has larger Poisson noise, hence 
unstable modes grow sooner. Note that the growth rate of stellar surface density in both 
models is in good agreement with the theoretical estimate for the fastest-growing mode, as 
represented by the long-dashed line in Figure HI The gaseous component, already driven 
nonlinear, has slightly faster growth. 
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Finally we report on the results of an important two-dimensional test. Evolution of 
linear waves in two-dimensional, shearing disks in principle would provide stringent tests 
for the current purposes, but we are unaware of test problems that have analytic solutions 
for comparison with numerical results. Swing amplification in a combined disk would also 
be a useful test, but it is difficult to accurately measure the amplification magnitude from 
numerical simulations. Propagating non-axisymmetric waves in a shearing disk also would 
be an excellent test problem for the gas component (cf. Paper I), but because of the La- 
grangian character it does not work as well for a stellar disk. We thus consider evolution 
of non-axisymmetric perturbations in non-shearing disks, for which equation (jSJ) still holds 
true if q — and k = (k 2 + k 2 ) 1 / 2 are used. We have run a set of numerical simulations 
for stable/unstable waves in various types of non-shearing disks, and confirmed that numer- 
ical results are consistent with the analytic predictions, verifying again the accuracy of our 
numerical code. 

Of particular interest is the effect of Poisson noise on gravitational instability in a 
stellar (or combined) system near the marginal state. Figure [5] plots the snapshots of surface 
density (in logarithmic scale) at t/t rb — 7.0 in a razor-thin, rigidly- rotating, star-only disk 
with q = 0, Q s = 1.2, and X = 2. Although the linear theory suggests that a disk with 
these parameters should be stable, the maximum surface density in a model with N pt \ = 10 4 
grows secularly, prompting the formation of 4 clumps, as the left panel of Figure [5] shows. 
The growth of perturbations in a model with N p t\ = 10 5 is slower, displaying over-density 
honeycomb structures (middle panel), while a model with N pt \ = 10 6 has a stable density 
field oscillating mildly (right panel). The mass- assignment scheme in our two-dimensional 
PM code yields Poisson noise at the level ST, s0 /'E s0 = 1.44(iV pt i/10 4 )~ 1//2 , as measured by the 
standard deviation £E s0 of the initial surface density S s0 - Some fraction of this noise may 
add to S s0 , decreasing the effective value of Q s . In cases when N pt \ is small and when Q s 
is slightly larger than unity, Q Sfi s can be smaller than the critical value and perturbations 
may grow artificiall;y0. This implies that Poisson noise can erroneously alter the fate of a 
numerically-modeled stellar system, especially when it is close to being marginal. 

This poi nt has sometimes been overlooked in previous iV-body simulations of disks. 



For example, iGriv et al.l (119991 ) simulated dynamics of non-shearing, local stellar disks using 
about 6,500 particles and found that disks with 1 < Q s0 £1.5 were unstable. They argued 
that the apparent d i screp ancy between the numerical and theoretical results might be partly 



because iLin fc Shul (1l966h 's WKB theory for density waves is incomplete and partly because 



2 If a half of 5S s o contributes directly to S s o, the effective Toomre parameter becomes Q s , e s = <3s(l + 
5S s0 /2E s0 ) _1 . When Q s = 1.2, Q SyCff = 0.70, 0.98, and 1.12 for iV pt] = 10 4 , 10 5 , and 10 6 , respectively, which 
is roughly consistent with the results of our numerical experiments. 
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their method for computing self-gravity is not very accurate. We suggest that large Poisson 
noise alone, simply associated with a small number of particles, can explain their numerical 
results. Adequate resolution is particularly an issue for global disk models, since the Poisson 
noise on a giv en spatial scale is proport i onal t o the linear size of th e disk. Some of the global 
simulations of lli. Mac Low, fc Klessenl koO^ \ and lLi et all (J20o3) (with 0.2 - 1.2 x 10 6 star 
particles for the whole disk), for example, may be affected by amplification of noise. 



5. Nonlinear Simulations 

To investigate nonlinear evolution of two-component differentially-rotating disks with 
a range of gas mass fractions, we have run a number of non-axisymmetric models. We fix 
q = 1, Q s — 2.1, X = 2, nH s /a S)X = 0.4, KH g /c g = 0.87, and vary Q g in the range of 
1 to 3. For these parameters, axisymmetric perturbations are stable. Solutions for self- 
consistent vertical equilibria show that H g decreases only by ~ 20% as Q g decreases from 



1.4 to 1.0 (IKim. Ostriker. fc Stond 120021 ) . so that we expect using a fixed value of H g does 



not significantly impact our quantitative results. 

We set up initially uniform gaseous and stellar disks, and apply density perturbations 
only to the gaseous component. Due to randomness in particle locations, the stellar part 
already has white-noise density perturbations with amplitudes of ~ 10% when N pt \ = 2 x 10 6 
which will be immediately transmitted into the gaseous disk. In order to keep the Poisson 
noise from dominating the evolution of the gaseous disk, we impose rather strong gaseous 
density perturbations using a Gaussian random field with a power spectrum (|£g ifc |) oc k~ 8 ^ 3 
for 1 < kL/(2ir) < 256 in Fourier space. This corresponds to a two-dimensional Kolmogorov 
spectrum if wave motions obey a sonic dispersion relation (Paper I). We measure the standard 
deviation of the gas surface density in real space and fix it to 10%. Although this level of 
density fluctuations is chosen as a compromise between maintaining the initial perturbations 
in the linear regime and reducing the effect of Poisson noise on the simulation outcomes, it 
may in fact well represent initial conditions for the highly-turbulent ISM observed in real 
disk galaxies. 



5.1. Threshold for Gravitational Runaway 

We begin with a high-resolution model with Q g = 1.0 and N p t\ = 2 x 10 6 . For this 
model, the combined gravity of gas and stars induces efficient swing amplification, leading 
to gravitational runaway. Because disk thickness dilutes self-gravity, each disk (i.e. purely 
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gaseous or purely stellar) in isolation would have remained stable both to axisymmetric 
and non-axisymmetric perturbations. Figure [6] plots as dashed lines the evolution of the 
maximum gaseous surface density and the spatially-averaged, stellar Toomre parameter Q s 
of this model (together with those of other models). Figure [7J displays snapshots of gas and 
stellar surface density in this model at three time epochs, while Figure [8] shows the initial 
modal growth of perturbations in the gas surface density over time. 

Initially, density perturbations (a superposition of modes with different wavenumbers) 
adjust by launching sound waves in the gas and redistributing stellar particles on their 
epicyclic orbits. During this relaxation phase, which lasts ~ 0.2t or b, coherent perturbations 
in stars and gas begin to develop through gravitational interactions. Due to the uniform 
background shear, the radial wavenumber k x (t) of perturbations increases linearly with time 
as 

k x (t) = k x (0) + qtlokyt, (9) 

where k y denotes the y- wavenumber and k x (0) is the initial x- wavenumber of perturbations. 
This kinematic shift of k x (t) causes wavefronts to swing from leading to trailing. Epicyclic 
motion of both stellar particles and gaseous fluid elements has the same rotational sense as 
the wavefronts, extending their exposure to self-gravity as they linger in wave crests. Pertur- 
bations k eep growing until the radial wavenumber becomes large, forming over-den se, trailing 



wavelets (IGoldreich fc Lynden-Belllll965bl ; I Julian fc Toomrdll966l ; lToomrelll98ll ; Paper I). 



The modal growth of density perturbations via swing amplification in the Q g = 1.0 
model is well illustrated in Figure M where we plot the evolution of power spectra of the 
gaseous surface density as functions of the normalized wavenumber n x (t) = Lk x (t) / (2ir) . 
Only the modes with k y = 2ir/L and |Ti a (0)| < 3 that have large amplitudes are shown. 
Note that a unity increment of n x (t) along each curve corresponds to a time elapse of 
At = (27r) _1 t or b- It is apparent that swing amplification is active only when waves are 
loosely wound with |na,(t)| ^3 — 4 above which sonic oscillations quench the perturbation 
growth. The n x (0) = 1 mode that happens to have the largest initial power grows and 
emerges first, as manifested by the trailing modes with n x (t/t or b = 0.3) ~ 3 shown in the 
top panels of Figure [71 However, it is the modes with n x (0) < that eventually dominate 
the initial swing amplification, since they have longer time for growth. 

As Figures [6] and [8] show, the initial swing amplification in the Q g = 1 model saturates 
at around t/t or b — 0.7. At this time, the gaseous surface density reaches £ 9;max /£ 90 = 18.7, 
while the stellar component reaches only S s max /E s0 = 2.5. The stellar component's density 
increases less due to its larger initial values of Q s and H s , and the steady heating associated 
with particle scattering off shearing wavelets. While the maximum gas surface density is 
rather high when the swing amplifier saturates, the gas filaments are very thin. Reduction of 
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self-gravity due to finite disk thickness prevents these filaments from immediately undergoing 
gravitational collapse. The filaments expand slightly, reducing their surface density. Some 
filaments have sufficient perturbed velocities to offset the tendency for wave rotation of the 
background shear. Consequently, the pitch angles of these high-density filaments vary little; 
they can be viewed as transient (swing-generated) spiral density waves in a two-component 
disk. 

Figure M plots sample cut profiles of various quantities along the x/L = —0.14 dotted 
lines in the snapshot shown in the middle row of Figure [7J The gas filament has a pitch 
angle % = 16° and is almost corotating with the center of the box; it coincides with a local 
stellar enhancement such that the composite may be considered a transient, local spiral 
arm segment. As the gas at —0.18 < x/L < —0.2 enters the arm, it is decelerated and 
compressed, analogous to a supersonic de Laval flow. This in turn increases the velocity 
paral l el to the arm due to the constraint of potential vorticity conservat ion (e.g., iHunter 
19641 iBalbus & Cowiel Il985l : iGammie! Il996l : iKim. Ostriker. fc Stone! bo02h . Note that the 
gravitational potential is dominated by the gaseous component , which tend s to symmetrize 
the density profile with respect to the peak (jLubow. Balbus. fc Cowiel Il986l ). The structure 
of this filament persists until t/t or ^ ~ 1.2 when collisions with other filaments make it highly 
self-gravitating. 

Paper I found that formation of bound clouds in swing-amplified disks occurs generically 
through one of the three secondary instabilities. In order of decreasing self-gravity, these are: 
(i) parallel fragmentation of filaments; (ii) gravitational collisions of shearing wavelets; and 
(iii) rejuvenated swing amplification preceded by nonlinear wave interactions. As explained 
above, high-density filaments in the Q g = 1.0 model do not immediately fragment due to 
thick-disk gravity. Only after several collisions with neighboring filaments is the density high 
enough to form gravitationally bound condensations. As Figure [7J shows, the Q g = 1.0 model 
forms 14 bound gaseous clouds, the average mass of which is ~ 6 x 1O 7 M , roughly 10 times 
larger than the Jeans mass Mj^d = c g/{G 2 ^ g o) in the corresponding gas-only, razor-thin 
disk. For comparison, in the Q g = 1 models of Paper I, the condensations that formed had 
M/Mj 2D ~ 0.5 — 1, suggesting that the effect of finite disk thickness (combined with an 
active stellar disk) is non-negligible in setting the cloud masses that form. 

We note that the bound condensations th at form in our models may represent H I su- 
percl ouds observed in many disk galaxies (e.g., lElmegreen fc Elmegreenlll983l ; iKnapen et al 



19931 ). In our models (isothermal, without feedback), these superclouds keep collapsing in 



a runaway fashion to the resolution limit. An artificial consequence of the runaway cloud 
collapse, at highly nonlinear stages, is the formation of loose stellar aggregates, as Figure 
[7J displays. Note that all the stellar clumps coincide with the massive gaseous clumps. In 
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realistic situations, however, gaseous superclouds would presumably fragment into lower- 
mass GMCs when appropriate physical processes such as internal turbulence are included. 
If M sc and R sc denote the mass and size of superclouds, respectively, stars would be drawn in 
strongly only if i? sc < GM SC /^ ~ 300 pc for M sc = 6 x 10 7 M Q and a s>x = 30 km s _1 . This 
corresponds to a supercloud surface density E sc > 200M Q pc -2 , comparable to the typical 
surface density of a GMC. As long as GMCs are destroyed by star formation before the whole 
supercloud reaches a density as large as that of an individual GMC (which is observationally 
true), stellar clumps would not be produced in reality. 

As shown in the right panel of Figure [6], the stellar velocity dispersion in this model 
continually increases as individual stars scatter off stellar and, more importantly, gaseous 



stella r disk by transient spiral dens i ty waves (e.g 
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20061 ). although the runaway fragmentation and collapse of gaseous fragments in our two- 
component models cause the stellar velocity dispersions to rise more rapidly than in star-only 
systems. 

Figure [6] shows that a model with Q g = 1.2 also becomes unstable, but with weaker 
gaseous gravity it takes longer to reach gravitational collapse. The first-generation filaments 
formed in this model do not fragment directly, and unlike in the Q g = 1 model collisions do 
not yield merger-induced fragmentation. Instead, the structures nonlinearly interact with 
each other and s upply fresh, small- \k x \ modes that undergo subsequent swing amplification 
(Paper I; see also IFuchsl 120051 for its stellar analog). Four successive stages of "rejuvenated" 
swing amplification (from interactions of filaments of order unity amplitudes) are required to 
drive the Q g = 1.2 model into eventual gravitational runaway. The average mass of gaseous 
clouds that form is again ~ 10Mj 2 d- 

When we further reduce the gas surface density so that Q g = 1.4, the initial swing 
amplification yields a moderately self-gravitating state that never forms bound clumps over 
the course of the simulation. As Figure [10] shows, the gaseous component in this model 
is dominated by shearing wavelets (and the related velocity field). Mild rejuvenated swing 
amplification allows the peak gaseous surface density and velocities to grow steadily, reaching 
the nonlinear regime. The stellar velocity dispersions also grow, but the net increase is less 
than 10% at t/t or b = 8 (see Fig. [6]); the stellar disk also remains virtually uniform with 
only very low fluctuations in density. We regard this model as "marginal" because it would 
almost certainly end up with bound condensations if evolved over a sufficiently long time. 
In other models with Q g > 1.7, growth is so weak that perturbed surface densities remain 
in the linear regime at saturation of the first swing amplification. Rejuvenation of swing is 
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almost absent in these models, so that the combined systems can be considered nonlinearly 
stable. We conclude, therefore, that the threshold for gravitational instability in combined 
gas-star disks has Q g near 1.4. 

5.2. Gravity-Driven Turbulence 

We now examine the properties of turbulence generated by swing amplification in two- 
component disks. We define 6E^ [n = | JS g (v — v ) 2 dxdy/J2„ , representing the specific 
kinetic energy associated with turbulent velocities in the the gasjj We have shown in §4.2l that 
large Poisson noise in a particle distribution effectively increases the mean surface density 
of stars in local regions, and can even incur spurious gravitational instability. Similarly, 
random particle noise can artificially - and possibly significantly - enhance the saturated 
state value of SE^. 

We want to measure SE^ in a way that minimizes the effect of particle noise. To 
this end, we take the following steps: (i) Run a two-component model in which both stellar 
and gaseous disks evolve under the total (gaseous plus stellar) self-gravity, and measure 
the resulting kinetic energy, which we denote by SE^. (ii) Run a control model with the 
same set of parameters as in step (i), but evolving the stellar component kinematically (i.e., 
omitting the V($ s + $ 9 ) gravity terms in eq. [I]), while evolving the gaseous component 
including both self- and stellar gravity. In runs of type (ii), no initial perturbation is applied 
to the gas, but the gas responds to the Poisson noise of the stars. The energy resulting 
from runs of type (ii) are denoted 6E^ . (iii) Calculate the gaseous velocity dispersion 
o v = (2(SE^) — 2(SE^)) 1 ^ 2 , where the brackets ( ) denote a temporal average over some 
time interval (typically over 1-6 orbits), (iv) Finally, vary iVp t i and repeat steps (i)-(iii) to 
check the convergence of a v with increasing N pt \. 

Figure [11] displays evolutionary histories of SE^ for Q g = 1.4 models that differ in 
the number of particles and the treatment of gravity for the stellar component. The thick 
curves correspond to runs following (i) above, while thin lines correspond to runs following 
(ii). In the cases with full gravity, the systems relax after the initial swing amplification, 
with SE^ decreasing temporarily as some of the gaseous kinetic energy is transferred via 
sound waves and/or weak shocks into thermal energy that dissipates under the isothermal 
prescription. The continuous input from stellar perturbations as well as nonlinear feedback 
from the shearing wavelets cause a resumed increase in SE^[. After t = 2 orbits, models 



3 Note that the total perturbed kinetic energy also contains other terms; SE^ is (half of) the mass- 
weighted mean-squared turbulent velocity dispersion. 
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with gravitating stars show a secular increase of 6E ^ . On the other hand, in models where 
the stars evolve passively, the initial swing is weak enough that SE^ simply saturates at a 
constant (low) value. Other models (not shown) with Q g > 2.0 show that both 5E^y and 
SE^ are more or less constant after t/t rb = 2, because swing amplification is very mild. 

Figure fl2Ta) plots the results for a v against iV p ti as filled circles for various Q g . In 
computing a v , we take time averages of 5-Ekin over t/t or b = 1 — 6 for all the models except the 
unstable model with Q g = 1.2. For Q g = 1.2, time averages are taken only over t/t OTh = 1—4 
since afterwards the velocity field is mainly a response to discrete high-density filaments. For 
comparison, Figure [T27 a) also plots values of a v = (25 E^) 1 / 2 as open circles for selected 
values of Q g . The figure shows that for Q g < 1.7, o v and a v are in close agreement. More 
generally, it is clear that there are progressively smaller differences between a v and a v as Q g 
decreases, indicating that the contribution from swing-amplified Poisson noise to a v becomes 
increasingly unimportant at lower Q g . 

Even after SE^ - embodying the kinetic energy directly induced in the gaseous disk 
by the stellar Poisson noise - has been subtracted, Figure [T2] shows that a v still depends 
on iVp t i. This is because the gas responds to swing- amplified perturbations in the stars that 
partly originate as noise in the initial stellar distribution. Note, however, that a v clearly 
converges with increasing N pt \, implying that this secondary noise effect becomes negligible. 
The convergence study suggests that o v at N pt \ = 2 x 10 6 provides a clean measure of 
the turbulent energy that is uncontaminated by effects of Poisson noise. Our results also 
imply that numerical simulations of marginally-unstable galactic disks are not reliable if 
insufficiently many particles are employed for the collisionless part. 

Figure fl27 b) plots o~ v for two-component disks with iVpti = 2 x 10 6 as well as the veloc- 
ity dispersion for gas-only disks. For stable systems with Q g £ 2.0, the stellar component 
enhances the amplitudes of the gaseous random velocities by about a factor of 2 compared 
to the gas-only counterpart, although velocity dispersions are small overall in this regime. 
When the disk is marginal or unstable (Q g = 1.2 or 1.4), a v becomes comparable to the 
effective speed of sound, and is ~ 4 — 5 times larger than the value obtained for a gas-only 
disk. This suggests not only that the stellar disk can significantly affect driving of gaseous 
turbulence, but also that self-gravity combined with galactic shear is capable of inducing 
sonic-level gaseous turbulence under realistic disk conditions. 

To characterize the turbulence generated by swing amplification, we plot in Figure [131 
Fourier power spectra of the compressive and shear components of gaseous velocities defined 
by 

-Pcomp = |k ■ 5v k \ 2 and P shC ar = |k x 5v k \ 2 , (10) 
respectively. Here, 5v k denotes the Fourier transform of the perturbed velocity, v — v . Data 
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at t/torb = 2.5 of a fully-gravitating two-component model with Q g = 1.4 and iV p ti = 2 x 10 6 
are used for Figure [131 We have also calculated the density power spectra and confirmed 
that they follow the compressive velocity power spectra very closely. 

The velocity power spectra ar e clearly anisotropic, as is typical of turbulence in systems 
with strong background shear (e.g.. lHawley. Gammie. fc Balbudll995l ; lKim. Ostriker. fc Stone 
20031 ). Comparisons of Figure [T3"b.fe reveal that the ratio of total power in the shearing to 
compressive parts is about 1.5. This value is q uite small compared to ~ 4 — 10 found in 



recent simulations of hypersonic turbulence (e.g.. iBoldyrev et all 120021 ; IVestuto et al.ll2003l ). 
presumably reflecting the facts that (i) our models are subsonic or weakly supersonic at 
best, so that compressions are not strongly dissipated, and (ii) driving by self-gravity fa- 
vors pumping of compressive motions. Although the power spectra vary too much at low 
k to yield clean p ower-law indices and are affected by numerical dissipation at high k (e.g., 
Kim fc Ryull2005l ). our results are consistent with P oc k~ 4,5 in the inertial range. The cuts 
along the k x -axis exhibit a break near k x L/ (2tt) ~ 8 — 10, below which the power flattens 
and above which P com p oc k~ A and P s hear oc k~ 6 . About 85% of the turbulent power is in large 
scales modes with A < L/10, suggesting that the small-scale turbulent velocity dispersion 
is a small fraction of the total; this proportion could, however, differ in three-dimensional 
models. The flattening of power spectra at small k x and steepening at large k x in our models 
makes sense, considering that in turbulence driven predominantly by swing amplification the 
background shear tends to increase the k x value of any structure over time. In quasi-steady 
state, nonlinear interactions of shearing wavelets evidently supply fresh power into small- \k x \ 
modes at just the rate needed to compensate for the kinematic shift of k x due to shear, such 
that the shape of the power spectrum is nearly flat. While turbulent energy cascades to small 
scales appear to shape power spectra along the fcy-direction, shear and nonlinear feedback 
clearly dominate energy flows in the fc^-direction. 



6. Summary and Discussion 
6.1. Summary 

In this paper we have analyzed the nonlinear dynamical evolution of two-component 
galactic disks. One primary goal has been to extend our previous work (Paper I), which 
modeled gas-only disks, to evaluate thresholds for runaway gravitational instability and 
formation of bound clouds. The other primary goal has been to quantify the generation 
of gaseous turbulence in cases where gravitational runaway does not occur. The gaseous 
medium is isothermal and evolved by a hydrodynamic technique as in Paper I, whereas the 
newly implemented stellar component is represented by an iV-body system evolved by a 
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particle-mesh method. The two components interact through their respective gravitational 
potentials. Our simulations are two-dimensional, but incorporate (important) finite disk 
thickness effects in an approximate fashion (see §2.11) . The local framework we employ 
incorporates galactic rotational shear, tidal gravity, and Coriolis forces. The present models 
do not include effects of magnetic fields or externally-driven stellar density waves. 

Our main results follows: 

1. Allowing for nonzero disk thickness significantly stabilizes axisymmetric modes in 
two-component disks (see §[3]). For solar neighborhood conditions with Q s = 2.1, the critical 
Toomre parameter of the gaseous component would be Q g , C rit = 1-27 if the both disks 
were razor-thin. Using more realistic scale heights H g = 170 pc and H s = 330 pc for the 
gaseous and stellar disks, respectively, the critical value is reduced to Q gtC rit = 0.67, well 
below the observed value. Two-component, thick disks with Q g = 1.4 are axisymmetrically 
stable unless Q s < 0.8, which is again inconsistent with observed galactic conditions. Star 
formatio n thresholds at Qn th ~ 1-4 observed in the o uter regions of external massive disk 



galaxies (IKennicuttl Il989l . iMartin fc Kennicuttl 120011 ) are therefore not a consequence of 



axisymmetric gravitational instability. 

2. Our two-dimensional simulations show that two-component disks undergo gravita- 
tional runway when Q g < Q 9jCr it (see §5 .11) . For stellar parameters representing the solar 
neighborhood (Q s = 2.1 and H s = 330 pc), this nonlinear threshold occurs at Q gtC rit ~ 1-4. 
Disks with Q g below this threshold experience nonaxisymmetric gravitational instability, 
forming bound condensations of mass ~ 6 x 10 M each. These condensation masses are 
10 times larger than those that form in razor-thin, gas-only disks because larger scales are 
required for gravitational instability in thicker disks. We note that the critical value of Q g for 
non-axisymmetric instability is larger than that for axisymmetric instability by more than a 
factor of 2. 

3. In addition to forming bound condensations, swing amplification is able to generate a 
significant level of gaseous turbulence (see §5.21) . The active stellar component is important to 
this, with velocity dispersions a factor 2 — 5 larger than in the equivalent gas-only disks. Our 
simulation results show that when Q g = 1.2 — 1.6, corresponding to disks near the stability 
threshold, the density-weighted velocity dispersions of the gas amount to a v ~ (0.3 — l)c s . 
This suggests that swing amplification can serve to tap rotational and gravitational energy 
in feeding random motions in the ISM. The turbulence in our models is anisotropic and has 
slightly more energy in the shearing than the compressive motions. Below k x L/ (2ir) £ 8 — 10, 
the power spectrum is relatively flat, while it steepens at larger k. Most (~ 85%) of the total 
power is contained in large scale modes with A < L/10, although in fully three dimensional 
models the relative power in small and large scales could potentially change. 



-21 - 



4. Poisson noise in the positions of randomly placed particles produces initial density 
fluctuation amplitudes « 1.4(A^ pt i/10 4 )~ 1//2 for two-dimensional particle-mesh simula- 

tions with grid resolution 256 2 . For near-marginal systems we show that these density varia- 
tions can lead to spurious local instabilities and artificial fragmentation ( §4.2p . Poisson noise 
can also lead to overestimates of turbulent velocity fluctuations produced by swing. Insuffi- 
cient particle numbers can therefore seriously compromise the results of iV-body disk simula- 
tions. Care must be taken to test numerical convergence. Here, we find that N pt \ = 1 — 2 x 10 6 
is needed to avoid contamination by Poisson noise. 



6.2. Discussion 



In this work, we have found that the threshold for formation of gravitationally bound gas 
clouds occurs at Q 9lCr it ~ 1.4. This result is in fact nearly identical to the nonlinear threshold 
value Q 9iC rit = 1-3 we found in Paper I (for unmagnetized models), due to compensation 
between two effects that were not accounted for in Paper I. The current models allow for a 
(stabilizing) nonzero disk thickness, and include a (destabilizing) active stellar component; 
although each effect by itself makes a significant difference in Q 9jC rit, in net they nearly 
cancel each other. While the present models do not include magnetic fields, we expect that 
inclusion of magnetic effects would not significantly change the results; in Paper I, there was 
only a 10% change in Q 9iCr i t even with equipartition-strength large-scale magnetic fields. 

An open question is how sensitive Q 9 , C rit is to the precise parameters characterizing the 
stellar disk. To define a manageable parameter space, we have fixed the value of Q s to 2.1, 
the (initial) ratio of the stellar velocity dispersion to the gas sound speed to 4.3, and the 
ratio of stellar to gas disk thicknesses to 1.9. The values chosen for these parameters are 
based on solar-neighborhood observations, and thus are representative of conditions in the 
middle of the galactic disk for a late-type spiral. Although the stellar component's param- 
eters surely vary from the deep interior to the outer reaches of a galactic disk, constraining 
the relevant quantities obseryationally is in fact quite chall enging. Extragalactic studies 
( Ivan der Kruit fc Searldll98ll . Il982l ; Ide Grijs fc Peletierlll997l ) are consistent with the stellar 
disk thicknesses being independent of radius at least for late-type spirals. In this case, a 
(self-gravitating) stellar disk with a constant ratio of vertical to radial velocity dispersion 

1 /2 

will have a x oc £ s q . For an exponential stellar surface density S s0 oc exp(—R/R d ) with 
scale length R d , this would imply that the variation of Q s over 1 < R/Rd < 4 is less than 
30%. Furthermore, Q s varies over only a small range from one galaxy to another; iBottema 
(119931) found that Q, 9 = 2 — 2.5 fits the data well (consistent with theoretical predictions of 



e.g. 



Sellwood fc Carlbergill984l ). Thus, although we obtained the result Q g , cr it ~ 1.4 for a 
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specific set of stellar parameters, the modest variation in these parameters among and within 
real galaxies suggests that there would be a correspondingly modest variation in Qg )Cr it- 

It is interesting to note that measured values of Q„ at the observed threshold rad ii 



for star formation are around Q 9t th ~ 1 — 2 (IKennicuttl 1 1983 ; iMartin fc Kennicuttl l200ll ). 



Unless the stellar properties at these threshold radii differ dramatically from those we have 
adopted, the results of our work, showing Q 9t th within the same range, quantitatively support 
the possibility that star formation thresholds and (nonaxisymmetric) gravitational instability 
thr esholds are likely to be o n e and the same. A similar conclusion was reached very recently 



by iLi. Mac Low, fc Klessenl (120051 ). who performed global SPH + iV-body simulations and 
found that a critical "combined star-gas" Toomre parameter of Q sg = 1.6 defines the limit 
between gravitationally-stable and unstable disks. Of course, turbulence driven by swing 
amplification could change the effective value of Q g if it contributes a stabilizing pressure 
in a manner similar to the microphysical random motions, and if the turbulent amplitudes 
are large. Since the threshold radii are by definition where disks become stable, however, 
then Q g > Q g ^it in those locations and there would not be a significant contribution from 
swing-driven turbulence. 

We note that many physical processes not included in the current models may poten- 
tially change Q g>C rit for gravitational runaway. These include tur bulence pervasive in the 



ISM, and thermal and other dynamical instabilities. For instance, iKim. Ostriker. fc Stone 



( 120031 ) found that nonlinear density fluctuations driven by magnetorotational instability in 
gas-only, three-dimensional disks enhance Q 9iCr it by about 50% relative to those in unmag- 
netized counterparts. It is especially important to explore the effects of turbulent motions 
and magnetic fields in multiphase systems, because the thermal velocity dispersion of cold 
gas is so low. The usual assumption is that an effective Q can be defined by appropriate 
weighting of various contributin g effective pressures and surface density components (cf. 



Wada. Meurer. fc Norman! 120021 ) ; an important direction for future research is to test this 
idea rigorouslyQ In future work, it will be interesting to test whether realistic turbulent, 
multi-phase models (in 3D or with thick-disk 2D gravity) show nonlinear threshold behavior 
for a suitably defined Q 9te s, and how a live stellar component quantitatively affects critical 
Q values. 

The fact that the stellar disk plays a significant role in destabilizing the gas disk may help 
explain the very low star formation rates in late-type, low surface brightness (LSB) galax- 



4 Single-phase models such as the present one are not useful for evaluating the idea or quantifying the 
value of Qg, e ff , since naive inclusion of the radial turbulent velocities that are present (< 0.5c g ) at both small 
and large scales increases Q g by < 10%. 
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ies. While LSB galax ies are often comparable in total gas mass to normal galaxies (e.g., 
Matthews et al.ll2005l ). most of them have gas surface densities below the threshold values 
^g,th corresponding to Q g = 1.4, such that low star formation rates and weak stellar disks 



are e asily explained (e.g., Ivan der Hulst et al.lll993l ; Ide Blok et al.lll996t lUson k Matthews 
20031 ) . However, some LSB galaxies with Q g < 1.4 still evidence little star formation 



([Pickering et al.lll99 7. 1999 :jO'Neil. Bothun. k Schombertl^OOOUO'Neil. Verheiien. k McGaugh 
2000l ; see also lElmegreenll2002l ) . In part, this is probably because Q s is so large that the stel- 
lar disk does little to encourage instability in the gas disk; we found that pure-gas disks are 
still quite stable when Q g = 1.2. In addition, without the vertical gr avity of a massive stellar 



disk, the gas disk thickness will increase by a factor of 1.5-2 (e.g. iKim. Ostriker. k Stone 
20021 ) compared to that in a normal galaxy. A larger value of H g dilutes gravity, which tends 
to further lower Q fl)Cr it. Together, these effects will reduce Q 9jCr it below 1 for LSB galaxies, 
making them more stable than previously thought. 



Galactic differential rotation represents a bountiful store of kinetic energy (e.g.. lvon Weizsacker 

195ll ). and we have shown that swing amplification in two-component disks is able to 
transform this energy into ISM turbulent motions with appreciable amplitudes. Other 
mechanisms that can channel she ared rotation into turbulence in the diffuse ISM include 
the magnetorotational instabilit y (jSellwood and Balbuslll999l ; IKim. Ostriker. k Stone! 120031: 



2004 



Piontek k Ostrikerl 12004 
Gomez k Cox! 



200 5]), and interactions with spiral arms ([Martos k Coxi Il998l ; 



Kim k Ostrikerl 120061 ; IKim. Kim, k Ostrikerl 120061 ) . Our current simu- 



lations show that provided that the disk is marginal or unstable, swing amplification is at 
least as efficient as the magnetorotational instability under conditions of comparable surface 
density. Two-phase self-gravitating hydrodynamic models dWada. Meurer. k Normanll2002r) 
and self-gravitating dissipative cloud- fluid models (IHuber k Pfennigerll200ll ) have similarly 
concluded that self-gravity coupled with galactic rotation can produce ISM turbulence that 
is sustained over many galactic rotations. 



Driving of turbulence at levels a v ^ 0.6c g requires a disk to be at least marginally unsta- 
ble. Because the outer, low-surface density regions of disks are quite stable, this mechanism 
cannot drive turbulence there@ But in the inner parts of disks, self-gravitational stirring in 
some regions (having slightly lower density) may go hand-in-hand with formation of massive, 
star- forming clouds in other regions (having slightly higher density). Because the fraction 
of gas converted to stars per cloud formation epoch is quite small ( 10%), the rate of 
reduction in the mean gas surface density - and hence rate of increase in Q g - is very low. 



5 The outer parts of disks, however, are also where the scale height flares; together with low H I surface 
density this implies low gas volume densities, which is most favorable for the magnetorotational instability 
to develop (Piontek & Ostriker 2006, in preparation). 
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As Q g increases over time, the timescale for cloud formation via gravitational instability 
also increases. The overall evolution is toward an asymptotic state in which Q g approaches 
Qg,crit- Galaxies may naturally "fine-tune" their parameters in such a way that self-gravity 
can help stabilize small scales - by generating turbulence - even as it forces large scales into 
collapse. 
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Fig. 1.— Marginal stability curves in the Q s -Q g plane for axisymmetric gravitational 
instability in two-component disks with R = c g /a x . The heavy curves correspond to thick 
disks with H g = 170 pc and H s = 330 pc, while the light curves are for razor-thin disks. The 
regions below each marginal curve represent conditions for which two-component disks are 
axisymmetrically stable. The solid circle at (Q s ,Q g ) = (2.1,1.4) marks solar neighborhood 
parameter values which are very stable when the disk is treated as thick, but would be barely 
stable if the disk were assumed razor-thin. 
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Fig. 2. — (Left) Critical Q g values and (right) the marginal perturbation wavelength A marg as 
functions of the gaseous scale height H g . Regions with Q g < Q gtCT i t or A > A marg correspond 
to a unstable parameter set. The solid lines correspond to the gas-only systems, while other 
lines are for combined disks with Q s = 2.1. The dots in both panels represent the solution 
for solar neighborhood parameters. 
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Fig. 3. — Oscillation frequencies of axisymmetric stellar density waves. The abscissa is the 
wavenumber k normalized by k SfiTit = k 2 /(27tGS s ). One-dimensional, razor-thin, star-only 
disks with Q s = 1.1 are adopted and iV p ti = 10 4 particles are used. Open circles and solid 
triangles are from the X = 2 and X = 1 models, respectively, both of which are in good 
agreement with the analytic solutions (lines) of equation ([8]). 
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Fig. 4. — Time evolution of maximum surface densities of gas and stars in a combined 
disk with q — 1, Q g — 1.2, Q s = 1.5, iJ g = 0.05c s //t, if s = 0.1cr x //c, and c^/o^ = 0.4. An 
adiabatic equation of state with 7 = 5/3 is assumed for the gas. The long-dashed line gives 
the theoretical growth rate of ~ 0.894fi 0) which is in good agreement with the results of 
numerical simulations. 
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Fig. 5. — Snapshots of surface density at t/t or b = 7.0 in non-shearing, star-only systems 
with Q s = 1.2, X = 2, H s = 0, but with differing particle numbers: (left) N ptl = 10 4 ; 
(middle) N pt \ = 10 5 ; (right) N pt \ = 10 6 . Gray-scale bars give log(S s /S s0 ). It is apparent 
that large Poisson noise associated with small iVp t i effectively increases the surface density 
locally, lowering Q s and thus destabilizing systems that would otherwise remain stable. 
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Fig. 6. — Evolution of (left) the maximum gaseous surface density and (right) the spatially- 
averaged stellar Toomre parameter Q s , for models with various Q g . All the models are two- 
dimensional and employ 2 x 10 6 stellar particles with <5 s ,init = 2.1. Models with Q g < 1.2 
experience gravitational runaway, while models with Q g > 1.7 are stable; the Q g = 1.4 model 
is marginal with a rapidly fluctuating density field. For unstable or marginal models, Q s 
increases considerably over time due to the gravitational interactions of stellar particles with 
both gaseous and stellar concentrations. 
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Fig. 7. — Snapshots of (left) gaseous and (right) stellar density in logarithmic color scale 
for the Q g = 1.0 model. Dotted lines in black in the middle row indicate the ^-position at 
which slices of various physical quantities are taken in Figure EE Strong filaments develop in 
the gas, and weaker filaments in the stars, due to swing amplification (t = 1.0 orbits). The 
gaseous filaments subsequently gravitationally fragment into clumps (t = 2.0 orbits). 
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Fig. 8. — Modal growth of power spectra |S 9jfc | 2 of the gaseous surface density against the 
normalized wavenumber n x (t) = Lk x (t) / (2n) in the Q g = 1.0 model for t/t orh £ 1.1. A few 
selected modes with k y = 2n/L and |ra x (0)| < 3 are shown. These loosely wound waves 
exhibit growth via swing amplification, then saturation. Note that n x (t) increases linearly 
with time, with At = An x (t)t or b/(2n). The t — point is the farthest left on each curve. 
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Fig. 9. — Profiles of gaseous surface density, gaseous velocities, and gravitational potentials 
across a filament at x/L = —0.14 when £/t or b = 1, for the Q g = 1.0 model. The profiles are 
quali tatively similar to those in gala ctic spiral shocks when both gas and stars are present 



(e.g. 



Lubow. Balbus. fc Cowidll986l ). 
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Fig. 10. — Snapshots of (left) gaseous and (right) stellar density in logarithmic color scale 
for the Q g = 1.4 model. Note that nonlinear shearing wavelets are prevalent in the gaseous 
disk, while density fluctuations in the stellar disk are very weak. 
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Fig. 11. — Evolutionary histories of the random kinetic energy of the gas, SE kin = | J S 9 (v — 
Vo) 2 dxdy/T, g0 , in two-component disks with Q g = 1.4 and Q s = 2.1. Three cases with 
different i\T pt i are shown. The thick curves are for the models where the both stellar and 
gaseous components are fully self-gravitating, while the thin lines correspond to the cases in 
which the combined gravity applies only to the gaseous disk. The stellar disk in the latter 
case evolves passively, but still imposes gravitational forcing from Poisson noise on the gas 
disk. 
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Fig. 12. — Amplitudes of mean fluctuating gaseous velocities a v in two-component models 
with Q s = 2.1, H g = 0.87c 9 /k,, and H s = 0Act x /k, to test numerical resolution effects. 
(a) Filled circles show a v as a function of stellar particle number 7V ptl after subtracting the 
random kinetic energy induced by Poisson noise in the stellar distribution, while open circles 
give a V) the gas velocity dispersion without this correction, (b) Filled circles give a v for 
JVpti = 2 x 10 6 , and filled triangles are for gas-only systems. See text for details. 
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Fig. 13. — Amplitudes of the power spectra along the k x - and A^-axes of (a) compressive 
and (b) shearing parts of the gaseous velocity in a two-component model with Q g — 1.4 at 
t/torb = 2.5. Solid line, dotted line, dashed line, and dot-dashed lines indicate the modes with 
k y L/(27i) = 0, 1, 2, 3 as a function of k x (heavy lines), and the modes with k x L/(2ir) = 0, 1, 
2, 3 as a function of k y (light lines), respectively. The velocity power is clearly anisotropic 
and the shearing part dominates at large scale. Both compressive and shearing parts along 
the k x -&xis show a break near k x L/(2n) ~ 8 — 10 below which the power flattens. 



